#include "Polynomicon.h"
#include "FourierTron.h"
#include <iostream>


using namespace std;
using namespace SeamStress;


Polynomicon::Polynomicon(FourierTron* f, vector<Seamstress*>* ss, unsigned long int nt, mpfr_prec_t p) : pins(this, ss), prec(p), retempval1(0,prec), retempval2(0,prec)
{
  if(nt==0){nthreads=1;}
  else{nthreads=nt;}
  if(ss==NULL){nthreads=1;}
  else{if(nthreads > ss->size()){nthreads = ss->size();}}
  
  fourier = f;
  
  tempval1.set(0,0,prec);
  tempval2.set(0,0,prec);
  
  thread_vec1.assign(nthreads, tempval1);
  rethread_vec1.assign(nthreads, nullptr);
  rethread_vec2.assign(nthreads, nullptr);
}


Polynomicon::~Polynomicon()
{
  
}


void Polynomicon::calculateFactorials(unsigned long int N)
{
  if(factorial.size() >= (N+1)){return;}
  factorial.clear();
  factorial_inv.clear();
  
  unsigned long int fsize = 1;
  unsigned long int log_fsize = 0;
  
  while(fsize<N){fsize <<= 1;log_fsize++;}
  
  retempval1.set(1);
  factorial.assign(fsize + 1, retempval1);
  factorial_inv.assign(fsize + 1, retempval1);
  
  unsigned long int thread_count = 1;
  unsigned long int log_thread_count = 0;
  while(thread_count<nthreads){thread_count <<= 1;log_thread_count++;}
  if(thread_count>nthreads){thread_count >>= 1;log_thread_count--;}
  while((log_thread_count<<1) > log_fsize){thread_count >>= 1;log_thread_count--;}
  while(thread_count > (fsize/thread_count)){thread_count >>= 1;log_thread_count--;}
  
  thread_tot = thread_count;
  thread_var1 = (fsize)/(thread_count);
  pins.sewStraight(&Polynomicon::factorialInit, thread_tot);
  
  thread_var1 = (1<<(log_fsize - log_thread_count));
  thread_var3 = thread_var1>>log_thread_count;
  for(unsigned long int j=1;j<thread_count;j++)
  {
    retempval1.set(factorial[j*thread_var1]);
    thread_var2 = j;
    pins.sewStraight(&Polynomicon::factorialRun, thread_tot);
  }
  
  
  pins.sewStraight(&Polynomicon::factorialInv, thread_tot);
}


void Polynomicon::factorialInit(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  unsigned long int start = 1 + w*thread_var1;
  unsigned long int finish = start + thread_var1 - 1;
  rethread_vec1[w].set(1);
  factorial[start].set(start);
  for(unsigned long int i=(start+1);i<finish;i++)
  {
    __builtin_prefetch(&(factorial[i+1]), 1, 3);
    rethread_vec1[w].set(i);
    mul(rethread_vec1[w], factorial[i-1], factorial[i]);
  }
  rethread_vec1[w].set(finish);
  mul(rethread_vec1[w], factorial[finish-1], factorial[finish]);
}


void Polynomicon::factorialRun(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  unsigned long int count = (thread_var2*thread_var1+1+w);
  unsigned long int max = thread_var3 - 1;
  
  __builtin_prefetch(&(factorial[count]), 1, 3);
  for(unsigned long int i=0;i<max;i++)
  {
    __builtin_prefetch(&(factorial[count+thread_tot]), 1, 3);
    factorial[count].mul(retempval1);
    count += thread_tot;
  }
  factorial[count].mul(retempval1);
}


void Polynomicon::factorialInv(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  unsigned long int i = w + 1;
  __builtin_prefetch(&(factorial[i]), 0, 3);
  __builtin_prefetch(&(factorial_inv[i]), 1, 3);
  unsigned long int max = factorial.size() - thread_tot;
  while(i < max)
  {
    __builtin_prefetch(&(factorial[i+thread_tot]), 0, 3);
    __builtin_prefetch(&(factorial_inv[i+thread_tot]), 1, 3);
    factorial_inv[i].div(factorial[i]);
    
    i += thread_tot;
  }
  factorial_inv[i].div(factorial[i]);
}





void Polynomicon::shiftPoly(vector<cxfloat>& poly, const cxfloat& Delta, bool set_scale, long int log_scale)
{
  if(poly.size() < 4)
  {
    unsigned int osize = poly.size();
    tempval1.set(0,0);
    poly.resize(4,tempval1);
    shiftPoly(poly, Delta);
    poly.resize(osize, tempval1);
    return;
  }
  
  unsigned long int thread_count = nthreads;
  if(thread_count > poly.size()){thread_count = poly.size();}
  
  if(set_scale==false)
  {
    p2_scale=0;
    unsigned long int M = 1;
    while(M < poly.size()){M<<=1;p2_scale++;}
    p2_scale--;p2_scale=-p2_scale;
  }
  else{p2_scale = log_scale;}
  tempvecpointer1 = &poly;
  thread_var1 = poly.size();
  thread_tot = thread_count;
  pins.sewStraight(&Polynomicon::scalePolyPower2, thread_tot);
  
  p2_scale=-p2_scale;
  cxfloat Shift(Delta);
  Shift.mul_power_2(p2_scale);
  shiftPolyBrute(poly, Shift);
  
  tempvecpointer1 = &poly;
  thread_var1 = poly.size();
  thread_tot = thread_count;
  pins.sewStraight(&Polynomicon::scalePolyPower2, thread_tot);
}


void Polynomicon::scalePolyPower2(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  unsigned long int max = thread_var1 - thread_tot;
  
  long int tscale;
  
  unsigned long int i=w;
  
  while(i < max)
  {
    tscale = i*p2_scale;
    (*tempvecpointer1)[i].mul_power_2(tscale);
    i += thread_tot;
  }
  tscale = i*p2_scale;
  (*tempvecpointer1)[i].mul_power_2(tscale);
}


void Polynomicon::shiftPolyBrute(vector<cxfloat>& poly, const cxfloat& Delta)
{
  calculateFactorials(poly.size() - 1);
  
  if(tempvec2.size() < poly.size()){tempval1.set(0,0);tempvec2.resize(poly.size(), tempval1);}
  if(tempvec3.size() < poly.size()){tempval1.set(0,0);tempvec3.resize(poly.size(), tempval1);}
  
  tempvecpointer1 = &poly;
  unsigned long int thread_count = nthreads;
  if(((poly.size())>>1) <= thread_count)
  {
    thread_count = ((poly.size())>>1);
    if(thread_count>1){thread_count--;}
    else if(thread_count==0){thread_count=1;}
  }
  
  if((poly.size() - thread_count - 1)<thread_count){thread_count = (poly.size() - thread_count - 1);}
  thread_tot = thread_count;
  thread_var1 = poly.size();
  pins.sewStraight(&Polynomicon::setUPoly, thread_tot);
  
  for(unsigned long int i=0;i<poly.size();i++){poly[i].set(1,0);}
  scalePoly(poly, Delta);
  tempvecpointer1 = &poly;
  if((poly.size() - thread_count - 1)<thread_count){thread_count = (poly.size() - thread_count - 1);}
  thread_tot = thread_count;
  thread_var1 = poly.size();
  pins.sewStraight(&Polynomicon::setVPoly, thread_tot);
  
  multiplyPoly(tempvec2, poly, tempvec3);
  
  thread_tot = thread_count;
  thread_var1 = poly.size();
  tempvecpointer1 = &poly;
  pins.sewStraight(&Polynomicon::finalizeShift, thread_tot);
}


void Polynomicon::setUPoly(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  unsigned long int i = w;
  
  __builtin_prefetch(&(factorial[i]), 0, 1);
  __builtin_prefetch(&((*tempvecpointer1)[i]), 0, 1);
  __builtin_prefetch(&(tempvec2[thread_var1 - 1 - i]), 1, 1);
  unsigned long int max = thread_var1 - thread_tot;
  while(i < max)
  {
    __builtin_prefetch(&(factorial[i+thread_tot]), 0, 1);
    __builtin_prefetch(&((*tempvecpointer1)[i+thread_tot]), 0, 1);
    __builtin_prefetch(&(tempvec2[thread_var1 - 1 - (i+thread_tot)]), 1, 1);
    mul(factorial[i], (*tempvecpointer1)[i], tempvec2[thread_var1 - 1 - i]);
    
    i += thread_tot;
  }
  mul(factorial[i], (*tempvecpointer1)[i], tempvec2[thread_var1 - 1 - i]);
}


void Polynomicon::setVPoly(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  unsigned long int i = w + thread_tot + 1;
  
  unsigned long int max = thread_var1 - thread_tot;
  
  i = w + 1;
  __builtin_prefetch(&(factorial_inv[i]), 0, 1);
  __builtin_prefetch(&((*tempvecpointer1)[i]), 1, 3);
  while(i < max)
  {
    __builtin_prefetch(&(factorial_inv[i+thread_tot]), 0, 1);
    __builtin_prefetch(&((*tempvecpointer1)[i+thread_tot]), 1, 3);
    (*tempvecpointer1)[i].mul(factorial_inv[i]);
    
    i += thread_tot;
  }
  (*tempvecpointer1)[i].mul(factorial_inv[i]);
}


void Polynomicon::finalizeShift(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  unsigned long int i = w;
  
  __builtin_prefetch(&(factorial_inv[i]), 0, 1);
  __builtin_prefetch(&(tempvec3[thread_var1 - 1 - i]), 0, 1);
  __builtin_prefetch(&((*tempvecpointer1)[i]), 1, 1);
  unsigned long int max = thread_var1 - thread_tot;
  while(i < max)
  {
    __builtin_prefetch(&(factorial_inv[i+thread_tot]), 0, 1);
    __builtin_prefetch(&(tempvec3[thread_var1 - 1 - (i+thread_tot)]), 0, 1);
    __builtin_prefetch(&((*tempvecpointer1)[i+thread_tot]), 1, 1);
    mul(factorial_inv[i], tempvec3[thread_var1 - 1 - i], (*tempvecpointer1)[i]);
    
    i += thread_tot;
  }
  mul(factorial_inv[i], tempvec3[thread_var1 - 1 - i], (*tempvecpointer1)[i]);
}




void Polynomicon::reciprocalPowerSeries(unsigned long int nterms, const vector<cxfloat>& poly, vector<cxfloat>& series)
{
  if(poly.size() < 2)
  {
    if(poly.size() == 1){series.resize(1, tempval1);series[0].set(poly[0]);series[0].reciprocal();}
    return;
  }
  
  vector<cxfloat> temp_series[2];
  vector<cxfloat> temp_poly;
  vector<cxfloat> remainder;
  vector<cxfloat> temp_vector;
  
  unsigned int which=0;
  
  temp_series[0].resize(2, tempval1);
  temp_series[0][0].set(poly[0]);temp_series[0][0].reciprocal();
  sqr(temp_series[0][0], temp_series[0][1]);temp_series[0][1].mul(poly[1]);temp_series[0][1].neg();
  
  unsigned long int cur_poly_size = temp_poly.size();
  temp_poly.resize(2, tempval1);
  for(unsigned long int i=cur_poly_size;i<temp_poly.size();i++){temp_poly[i].set(poly[i]);}
  
  unsigned long int n = 1;
  
  while(temp_series[which].size() < nterms)
  {
    unsigned long int cur_poly_size = temp_poly.size();
    if(((n<<1) + 1) <= poly.size()){temp_poly.resize((n<<1) + 1, tempval1);}
    else{temp_poly.resize(poly.size(), tempval1);}
    for(unsigned long int i=cur_poly_size;i<temp_poly.size();i++){temp_poly[i].set(poly[i]);}
    if(((n<<1) + 1) > poly.size()){tempval1.set(0,0);temp_poly.resize((n<<1) + 1, tempval1);}
    
    
    multiplyPoly(temp_poly, temp_series[which], temp_vector);
    remainder.clear();
    for(unsigned long int i=(n+1);i<((n<<1)+1);i++){remainder.push_back(temp_vector[i]);}
    
    tempval1.set(0,0);
    temp_vector.clear();
    temp_vector.resize((n<<1) + 1, tempval1);
    temp_vector[0].set(1,0);
    for(unsigned long int i=0;i<remainder.size();i++){temp_vector[i + n + 1].set(remainder[i]);temp_vector[i + n + 1].neg();}
    
    multiplyPoly(temp_series[which], temp_vector, temp_series[(which+1)%2]);
    tempval1.set(0,0);
    temp_series[(which+1)%2].resize((n<<1) + 1, tempval1);
    
    which = (which+1)%2;
    
    n<<=1;
  }
  
  series.resize(nterms, tempval1);
  for(unsigned long int i=0;i<nterms;i++){series[i].set(temp_series[which][i]);}
}


void Polynomicon::reciprocalPowerSeries(unsigned long int nterms, const refloat_vec& poly, refloat_vec& series)
{
  if(poly.size() < 2)
  {
    if(poly.size() == 1){series.resize(1, retempval1);series[0].set(poly[0]);series[0].reciprocal();}
    return;
  }
  
  refloat_vec temp_series[2];
  refloat_vec temp_poly;
  refloat_vec remainder;
  refloat_vec temp_vector;
  
  unsigned int which=0;
  
  temp_series[0].resize(2, retempval1);
  temp_series[0][0].set(poly[0]);temp_series[0][0].reciprocal();
  sqr(temp_series[0][0], temp_series[0][1]);temp_series[0][1].mul(poly[1]);temp_series[0][1].neg();
  
  unsigned long int cur_poly_size = temp_poly.size();
  temp_poly.resize(2, retempval1);
  for(unsigned long int i=cur_poly_size;i<temp_poly.size();i++){temp_poly[i].set(poly[i]);}
  
  unsigned long int n = 1;
  
  while(temp_series[which].size() < nterms)
  {
    unsigned long int cur_poly_size = temp_poly.size();
    if(((n<<1) + 1) <= poly.size()){temp_poly.resize((n<<1) + 1, retempval1);}
    else{temp_poly.resize(poly.size(), retempval1);}
    for(unsigned long int i=cur_poly_size;i<temp_poly.size();i++){temp_poly[i].set(poly[i]);}
    if(((n<<1) + 1) > poly.size()){retempval1.set(0);temp_poly.resize((n<<1) + 1, retempval1);}
    
    
    multiplyPoly(temp_poly, temp_series[which], temp_vector);
    remainder.clear();
    for(unsigned long int i=(n+1);i<((n<<1)+1);i++){remainder.push_back(temp_vector[i]);}
    
    retempval1.set(0);
    temp_vector.clear();
    temp_vector.resize((n<<1) + 1, retempval1);
    temp_vector[0].set(1);
    for(unsigned long int i=0;i<remainder.size();i++){temp_vector[i + n + 1].set(remainder[i]);temp_vector[i + n + 1].neg();}
    
    multiplyPoly(temp_series[which], temp_vector, temp_series[(which+1)%2]);
    retempval1.set(0);
    temp_series[(which+1)%2].resize((n<<1) + 1, retempval1);
    
    which = (which+1)%2;
    
    n<<=1;
  }
  
  series.resize(nterms, retempval1);
  for(unsigned long int i=0;i<nterms;i++){series[i].set(temp_series[which][i]);}
}









void Polynomicon::shiftPoly(refloat_vec& poly, const refloat& Delta, bool set_scale, long int log_scale)
{
  if(poly.size() < 4)
  {
    unsigned int osize = poly.size();
    retempval1.set(0);
    poly.resize(4,retempval1);
    shiftPoly(poly, Delta);
    poly.resize(osize, retempval1);
    return;
  }
  unsigned long int thread_count = nthreads;
  if(thread_count > poly.size()){thread_count = poly.size();}
  
  if(set_scale==false)
  {
    p2_scale=0;
    unsigned long int M = 1;
    while(M < poly.size()){M<<=1;p2_scale++;}
    p2_scale--;p2_scale=-p2_scale;
  }
  else{p2_scale = log_scale;}
  retempvecpointer1 = &poly;
  thread_var1 = poly.size();
  thread_tot = thread_count;
  pins.sewStraight(&Polynomicon::scalePolyPower2Real, thread_tot);
  p2_scale=-p2_scale;
  refloat Shift(Delta);
  Shift.mul_power_2(p2_scale);
  shiftPolyBrute(poly, Shift);
  retempvecpointer1 = &poly;
  thread_var1 = poly.size();
  thread_tot = thread_count;
  pins.sewStraight(&Polynomicon::scalePolyPower2Real, thread_tot);
}


void Polynomicon::scalePolyPower2Real(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  unsigned long int max = thread_var1 - thread_tot;
  
  long int tscale;
  
  unsigned long int i=w;
  
  while(i < max)
  {
    tscale = i*p2_scale;
    (*retempvecpointer1)[i].mul_power_2(tscale);
    i += thread_tot;
  }
  tscale = i*p2_scale;
  (*retempvecpointer1)[i].mul_power_2(tscale);
}






void Polynomicon::shiftPolyBrute(refloat_vec& poly, const refloat& Delta)
{
  calculateFactorials(poly.size() - 1);
  if(retempvec2.size() < poly.size()){retempval1.set(0);retempvec2.resize(poly.size(), retempval1);}
  if(retempvec3.size() < poly.size()){retempval1.set(0);retempvec3.resize(poly.size(), retempval1);}
  
  retempvecpointer1 = &poly;
  unsigned long int thread_count = nthreads;
  if(((poly.size())>>1) <= thread_count)
  {
    thread_count = ((poly.size())>>1);
    if(thread_count>1){thread_count--;}
    else if(thread_count==0){thread_count=1;}
  }
  if((poly.size() - thread_count - 1)<thread_count){thread_count = (poly.size() - thread_count - 1);}
  thread_tot = thread_count;
  thread_var1 = poly.size();
  pins.sewStraight(&Polynomicon::setUPolyReal, thread_tot);
  
  for(unsigned long int i=0;i<poly.size();i++){poly[i].set(1);}
  scalePoly(poly, Delta);
  retempvecpointer1 = &poly;
  if((poly.size() - thread_count - 1)<thread_count){thread_count = (poly.size() - thread_count - 1);}
  thread_tot = thread_count;
  thread_var1 = poly.size();
  pins.sewStraight(&Polynomicon::setVPolyReal, thread_tot);
  
  multiplyPoly(retempvec2, poly, retempvec3);
  
  thread_tot = thread_count;
  thread_var1 = poly.size();
  retempvecpointer1 = &poly;
  pins.sewStraight(&Polynomicon::finalizeShiftReal, thread_tot);
}


void Polynomicon::setUPolyReal(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  unsigned long int i = w;
  
  __builtin_prefetch(&(factorial[i]), 0, 1);
  __builtin_prefetch(&((*retempvecpointer1)[i]), 0, 1);
  __builtin_prefetch(&(retempvec2[thread_var1 - 1 - i]), 1, 1);
  unsigned long int max = thread_var1 - thread_tot;
  while(i < max)
  {
    __builtin_prefetch(&(factorial[i+thread_tot]), 0, 1);
    __builtin_prefetch(&((*retempvecpointer1)[i+thread_tot]), 0, 1);
    __builtin_prefetch(&(retempvec2[thread_var1 - 1 - (i+thread_tot)]), 1, 1);
    mul(factorial[i], (*retempvecpointer1)[i], retempvec2[thread_var1 - 1 - i]);
    
    i += thread_tot;
  }
  mul(factorial[i], (*retempvecpointer1)[i], retempvec2[thread_var1 - 1 - i]);
}


void Polynomicon::setVPolyReal(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  unsigned long int i = w + thread_tot + 1;
  
  unsigned long int max = thread_var1 - thread_tot;
  
  i = w + 1;
  __builtin_prefetch(&(factorial_inv[i]), 0, 1);
  __builtin_prefetch(&((*retempvecpointer1)[i]), 1, 3);
  while(i < max)
  {
    __builtin_prefetch(&(factorial_inv[i+thread_tot]), 0, 1);
    __builtin_prefetch(&((*retempvecpointer1)[i+thread_tot]), 1, 3);
    (*retempvecpointer1)[i].mul(factorial_inv[i]);
    
    i += thread_tot;
  }
  (*retempvecpointer1)[i].mul(factorial_inv[i]);
}


void Polynomicon::finalizeShiftReal(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  unsigned long int i = w;
  
  __builtin_prefetch(&(factorial_inv[i]), 0, 1);
  __builtin_prefetch(&(retempvec3[thread_var1 - 1 - i]), 0, 1);
  __builtin_prefetch(&((*retempvecpointer1)[i]), 1, 1);
  unsigned long int max = thread_var1 - thread_tot;
  while(i < max)
  {
    __builtin_prefetch(&(factorial_inv[i+thread_tot]), 0, 1);
    __builtin_prefetch(&(retempvec3[thread_var1 - 1 - (i+thread_tot)]), 0, 1);
    __builtin_prefetch(&((*retempvecpointer1)[i+thread_tot]), 1, 1);
    mul(factorial_inv[i], retempvec3[thread_var1 - 1 - i], (*retempvecpointer1)[i]);
    
    i += thread_tot;
  }
  mul(factorial_inv[i], retempvec3[thread_var1 - 1 - i], (*retempvecpointer1)[i]);
}


void Polynomicon::addPoly(const vector<cxfloat>& to_add, vector<cxfloat>& added_to)
{
  tempval1.set(0,0);
  if(to_add.size() > added_to.size()){added_to.resize(to_add.size(), tempval1);}
  
  thread_tot = nthreads;
  if(thread_tot > to_add.size()){thread_tot = to_add.size();}
  
  tempvecpointer2 = &to_add;
  tempvecpointer1 = &added_to;
  
  thread_var1 = to_add.size();
  pins.sewStraight(&Polynomicon::addPolyThread1, thread_tot);
}


void Polynomicon::addPolyThread1(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  unsigned long int part = thread_var1/thread_tot;
  unsigned long int rem = thread_var1 - part*thread_tot;
  unsigned long int start = 0;
  unsigned long int end = 0;
  if(w<rem)
  {
    start = (part+1)*w;
    end = start + part;
  }
  else
  {
    start = (part+1)*rem;
    start += (part*(w-rem));
    end = start + part - 1;
  }
  
  for(unsigned long int i=start;i<=end;i++)
  {
    (*tempvecpointer1)[i].add((*tempvecpointer2)[i]);
  }
}


void Polynomicon::addPoly(const refloat_vec& to_add, refloat_vec& added_to)
{
  retempval1.set(0);
  if(to_add.size() > added_to.size()){added_to.resize(to_add.size(), retempval1);}
  
  thread_tot = nthreads;
  if(thread_tot > to_add.size()){thread_tot = to_add.size();}
  
  retempvecpointer2 = &to_add;
  retempvecpointer1 = &added_to;
  
  thread_var1 = to_add.size();
  pins.sewStraight(&Polynomicon::addPolyThread1Real, thread_tot);
}


void Polynomicon::addPolyThread1Real(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  unsigned long int part = thread_var1/thread_tot;
  unsigned long int rem = thread_var1 - part*thread_tot;
  unsigned long int start = 0;
  unsigned long int end = 0;
  if(w<rem)
  {
    start = (part+1)*w;
    end = start + part;
  }
  else
  {
    start = (part+1)*rem;
    start += (part*(w-rem));
    end = start + part - 1;
  }
  
  for(unsigned long int i=start;i<=end;i++)
  {
    (*retempvecpointer1)[i].add((*retempvecpointer2)[i]);
  }
}


void Polynomicon::subPoly(const vector<cxfloat>& to_subtract, vector<cxfloat>& subtracted_from)
{
  tempval1.set(0,0);
  if(to_subtract.size() > subtracted_from.size()){subtracted_from.resize(to_subtract.size(), tempval1);}
  
  thread_tot = nthreads;
  if(thread_tot > to_subtract.size()){thread_tot = to_subtract.size();}
  
  tempvecpointer2 = &to_subtract;
  tempvecpointer1 = &subtracted_from;
  
  thread_var1 = to_subtract.size();
  pins.sewStraight(&Polynomicon::addPolyThread1, thread_tot);
}


void Polynomicon::subPolyThread1(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  unsigned long int part = thread_var1/thread_tot;
  unsigned long int rem = thread_var1 - part*thread_tot;
  unsigned long int start = 0;
  unsigned long int end = 0;
  if(w<rem)
  {
    start = (part+1)*w;
    end = start + part;
  }
  else
  {
    start = (part+1)*rem;
    start += (part*(w-rem));
    end = start + part - 1;
  }
  
  for(unsigned long int i=start;i<=end;i++)
  {
    (*tempvecpointer1)[i].sub((*tempvecpointer2)[i]);
  }
}









